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We report a real-space renormalization group (RSRG) algorithm, which is formulated through 
Baxter's corner transfer matrix (CTM), for two-dimensional (d — 2) classical lattice models. 
The new method performs the renormalization group transformation according to White's 
density matrix algorithm, so that variational free energies are minimized within a restricted 
degree of freedom m. As a consequence of the renormalization, spin variables on each corner 
of CTM are replaced by a m-state block spin variable. It is shown that the thermodynamic 
functions and critical exponents of the q = 2, 3 Potts models can be precisely evaluated by use 
of the renormalization group method. 
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§1. Introduction 

The density matrix renormalization group (DMRG) established by WhiteSi has greatly en- 
hanced the applicability of the numerical real-space renormalization group (RSRG)iBi to one- 
dimensional (d = 1) quantum systems. After the applications to the 5 = 1 Heisenberg spin 
chainJlUB DMRG has been applied to a number of d = 1 quantum systems, such as Heisenberg 
ladder,!'!' 00 bond-alternating spin chain, 0'Ei§ strongly correlated electron systemJEHil'El'il) 
etc. On one hand improvements upon the numerical algorithm of DMRG have been proposed in 
order to analyze impurity systems JlBIl^ random systems, Bethe lattice systems, EH correlated 
electron system defined in the momentum space,© spin chains under finite temperature,© and so 
on. 

Recently White has analyzed the ground state wave function of d = 2 quantum spin systems, using 
a numerically accelerated finite system DMRG algorithm.© At first, White's acceleration technique 
was only applicable to the finite system DMRG algorithm. Soon after the authors extended the 
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acceleration technique to infinite system DMRG algorithm.© Further numerical acceleration in 
DMRG has been reported for models that have (quantum-) group symmetry. 

The DMRG picks up relevant correlations between the local system (=block) and the rest of the 
system (=reservoir). Since irrelevant spin (or particle) fluctuations are projected out, the DMRG 
method intrinsically has a variational property. Indeed, the ground state wave function obtained 
by DMRG is a good variational wave function, which is written as a product of tensors. BBE3 
In the thermodynamic limit the tensors become position independent, where the DMRG coincides 
with Ostlund's variational method. Martin-Delgado, Rodriguez-Laguna and Sierra have re- 
formulated the variational relation in DMRG using projection operators, and have proposed new 
RSRG algorithms. SB© 

Since d = 1 quantum systems are naturally related to d = 2 classical systems, it is possible to 
apply DMRG algorithm to the latter© The largest eigenvalue of the row-to-row transfer matrix 
is primarily important for the analysis of the d = 2 classical models. The DMRG applied to a 
d = 2 classical model evaluates the lower bound of the largest eigenvalue, using a variational state 
vector written in a product of tensors. It has been shown that the thermodynamic functions can 
be obtained precisely by DMRG in off critical regions. However, the numerical convergence in 
free energy becomes slow near the critical temperature, and therefore extensive computations are 
necessary at criticality. The reason of this slowing down is that the maximum eigenvalue of the 
row-to-row transfer matrix is nearly degenerate in the critical region. Such a degeneracy spoils 
the numerical efficiency of the Lanczos diagonalization, that is the most time consuming part in 
DMRG. 

Baxter's method of corner transfer matrixIzHill © (CTM), that was formulated in 1968 as an 
extension of the Kramers- Wannier approximation, 

0-0) 

is another variational method for d = 2 
classical lattice systems. The method gives approximate free energy per site in the thermodynamic 
limit. Baxter's method seems to have no relation with DMRG, but actually both of them are 
deeply connected; they are both iterative renormalization group method, and have the same fixed 
point in the thermodynamic limit. Baxter's method can be used as a numerical met hod ^3 EEH) 
and it runs faster than DMRG at criticality. This is because the largest eigenvalue of CTM is not 
degenerate even at the critical point. 

In this paper we introduce the advantage of Baxter's method into the numerical algorithm of 
DMRG for d = 2 classical system. We express the density matrix as a product of four CTMs. For the 
brevity, we call the improved renormalization group method 'corner transfer matrix renormalization 
group' (CTMRG) in the following© Apart from Baxter's method, the purpose of CTMRG is to 
obtain variational free energies of finite size systems, up to a certain system size. 

In the next section we review the construction of CTM. In order to simplify the discussion, we 
consider a (/-state Potts model on a decorated lattice, because it is easy to define CTM on the 
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lattice. In §3 we explain the variational relation in DMRG using projection operators. We then 
introduce CTM into the formulation of DMRG, and make up the numerical algorithm of CTMRG 
in §4. In §5 we apply CTMRG to q = 2, 3 Potts models. It is verified via the finite size scaling 
analysis^'© that CTMRG gives correct critical exponents. Conclusions are summarized in §6. 

§2. Corner Transfer Matrix 

We consider a g-state Potts mode]El§) on a decorated square lattice, whose geometry is shown 
in Fig.l. The white marks represent g-state spins {s} on vertices, and the black ones represent 
another set of g-state spins {a} that are in between s-spins. We refer to the model as 'decorated 
Potts model' in the following. 



Fig. 1. The g-state Potts model on a decorated lattice. The black and white marks represent q-state spin variables 
{a} and {s}, respectively. 



We consider a square shaped finite size system of linear dimension L; the case L = 3 is shown in 
Fig.l. The partition function is 

M W { (ij) J 
where K* is the interaction parameter, (ij) specifies the neighboring as spin pairs, and 5(ai,Sj) 
is equal to unity if <7j = sj and zero otherwise. When q = 2 the model coincides with the super 
exchange Ising model by Fisher.il' Since the decorated lattice is bipartite, we can explicitly take the 
configuration sum over u-spins, leaving that for s-spins. The partition function after the summation 
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Z L =(q-2 + 2e K ') M ^MKY J Ksi,s j )\, (2.2) 

w I m J 

where M = 2L(L — 1) is the number of c-spins in the system of the linear dimension L, and K is 
the effective interaction parameter between s-spins; K is determined through the duality relation 
for the d = 1 Potts mo del© 

e K = 1 1 + e (2 .3) 
q - 2 + 2e^ V ; 

Equation (2.2) shows that Z L is proportional to the partition function of the Potts model on the 

simple square lattice. 

In the same manner, we can decimate s-spins and express Z L in eq.(2.1) as a partition function 
of a symmetric vertex model. The vertex weight, that represents four cr-spin interaction around a 
vertex, is expressed as 

W abcd = J2 ew^K*(5 sa + 5 sb + 5 SC + 5 sd ) j , (2.4) 

where s represents an s-spin on the vertex, and a,b,c,d represent the <r-spins around s. The weight 
W a hcd is invariant under the permutations of indices, since our model is isotropic. (See eq.(2.1).) 
We refer to the vertex model as 'g 4 -vertex model', since the indices a, 6, c, and d run from 1 to q, 
and since all of the weights are nonzero. 

In order to express Z L correctly, we have to define two additional Boltzmann weights at the 
boundary. One is 

P abc = ]Texp K*(6 sa + 5 sb + 5 SC ) \ (2.5) 



that expresses the three-spin interaction at the side of the square system, and the other is 

C ab = j2exp^K*(5 sa + 5 sb )^ (2.6) 

that represents two-spin interaction at the corner. The partition function Z L in eq.(2.1) is expressed 
as a tensor product of these Boltzmann weights; for example, the partition function when L = 3 is 

^ 3 = X] W kheb P abc C cd P de fC 'fgPg^C^Pj^Cia, (2.7) 

ab...l 

where the arrangement of spin indices is shown in Fig. 2. 

In order to generalize the expression of Z 3 in eq.(2.7) to larger systems, we introduce the half-row 
transfer matrix© (HRTM) and CTM. The HRTM is the left (or the right) half of the row-to-row 
transfer matrix, which is a generalization of the boundary weight in eq.(2.5). The HRTM of length 
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Fig. 2. The symmetric g 4 -vertex model that corresponds to the decorated q-state Potts model in Fig.l. 



N can be defined by the recursion relation 



P. 



>N _ 
abc — 



pN-1 
a N d cjv b r a'dC'i 



(2.8) 



d=i 



where the vector index 



a = (a±, a,2, ■ ■ ■ , ajv-i, ajv) 



(2.9) 



represents a group of (/-state spin indices on the half-row of length N, and a' = (a±, . . . , a^-i) is 
included in a = (o', ajv); the same for c = (c', cjv). The initial value Pa bc , where a = (a\) and 
c = (ci), is given by the boundary weight in eq.(2.5). As an example, we show the case iV = 3 in 
Fig.3. We often abbreviate the vector indices of P^bc an d write the HRTM as P b N ; in such a case 
we think of P b N as a ^^-dimensional matrix (P b N )ac- 

The CTM is a generalization of the boundary weight C a b in eq.(2.6). We define CTM using the 
recursion relation 



where we have used the index rule a = (a', a^r) in eq.(2.9); the same for b, c and d. The initial 
value C^, where a = (a\) and b = (b\), is given by the corner weight in eq.(2.6). Figure 4 shows 
the CTM when N = 3. The factor WP^^Pf' 1 inside the parenthesis in eq.(2.10) plays the role 
of a transfer matrix that extends the area of CTM. 




(2.10) 
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Fig. 3. The half-row transfer matrix of length 3 (= P 3 ), and its components P 2 and W. 



It should be noted that our definition of CTM in eq.(2.10) is different from the conventional 
one; C N in eq.(2.10) corresponds to a square quadrant of the whole system, while Baxter's CTM 



HRTMs are necessary when we increase the size of CTM. (See eq.(2.10)) 

We have so far chosen free boundary conditions. Since the boundary weights P a bc an d C a b in 
eqs.(2.5) and (2.6) determine the boundary conditions, we can impose fixed boundary conditions 
by modifying the definition of P a b c and C a b- For example, if we impose Sj = 1 for all i-site at the 
boundary, the boundary weights become 



corresponds to a triangular region of a square lattice system.EP Because of this difference, two 



P 





according to the fixed boundary condition. 

Now we can express the partition function Z L 
systems, the partition function is written as 



in terms of P and C . For even-size (L = 2N) 





abed 
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Fig. 4. Corner Transfer Matrix C 3 , and its components P 2 , C 2 , and W. The shaded region corresponds to C 2 . 



For odd-size (L = 2N + 1) cases, we generalize eq.(2.7) and express Z 2N+1 as 

Z 2N+1 = £ W kheb Tr(p k N C N P h N C N P e N C N P b N C N ) (2.13) 

kheb 

where we have regarded P N as a matrix. 

In addition to the partition function, we can obtain thermal averages of the spin polarization or 
spin correlation functions in the same way. For example, the s-spin at the center of the odd-size 
(L = 2N + 1) cluster takes the direction '1' with the probability 

v Trf pN /~<N pN /~iN pN /~tN pN /~<N\ 
v l^abcd^abcd^Ya U r b ° r c U r d ° ) 

(Sis) = 7 r> (2-14) 

E abcd W abcd Tr(p?C"P b "C"P?C"P?C") 

where X is a new vertex weight 

X abcd = exp jiT(5i a + S lb + 6 lc + (2.15) 

that counts the case s = 1. 

The corner transfer matrix is g^-dimensional, where the dimension increases rapidly with 
N. The fact prevents us from exact numerical calculation of Z L by way of eq.(2.12) and (2.13). 
For example when q = 2, the upper limit of N is about 13. The restriction for N is more sever for 
larger q cases. This is the main reason that we employ the renormalization group (RG) method. 
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§3. Minimum Free Energy Principle 

The background of DMRG is the minimum free energy (or the maximum partition function) 
principle which is represented by the density matrix. Generally speaking, a density matrix p in 
statistical mechanics is a matrix whose trace coincides with the partition function of the system. 
The product of four CTMs 

P ™=(C N ) A (3.1) 

is a kind of density matrix, whose trace Tr p 2N is the partition function Z in eq.(2.12). Figure 
5 shows the construction of p 2N when N = 3; the element of the density matrix p 2 ^ represents 
the Boltzmann weight of the system that has a cut (or gap), where a and b represent the spin 
configurations on each side of the cut. 



Fig. 5. Density matrix p 2N for N = 3. The shaded region represents the corner transfer matrix C 3 . 



Equation (3.1) shows that both C N and p 2N have the common eigenvectors Ri (1 < i < q N ) 
that satisfy the eigenvalue equations 

C N Ri = Ri LOi 

p 2N R i = R i ut, (3.2) 

where u>i and uf is the i-th eigenvalue of C N and p 2N , respectively. We write down the element of 
Ri as R a j in the following. According to the symmetry of the g 4 -vertex model, both C N and p 2N 
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are symmetric. Therefore the square matrix R = (R\, R2, • • •) satisfies the orthogonal relation 

[R t R) =(Ri,Rj) = Sij 
(R RT ) ab = E R aiR bi = d ab- (3-3) 

i 

We assume the decreasing order w 1 4 > w 2 4 > • ■ for the eigenvalues, where all of the Wj 4 are 
nonnegative; as far as our q -vertex model is concerned, oj 2 is also real and nonnegative. A trivial 
inequality 

m 

Z2N '^X>$ 4 (3.4) 
£=1 

holds for arbitrary integer m; the l.h.s. coincides with the r.h.s. when m = q N , or when uji = for 
all i > m. (Hereafter we use greek letters for indices that run from 1 to m.) It has been known 
that the dumping of ujf with respect to i is fairly rapid. Therefore the r.h.s. of eq.(3.4) is 

a good approximation for Z 2N for sufficiently large m. Typically m is of the order of hundreds in 
realistic numerical calculations. 
Let us rewrite eq.(3.4) into the matrix formula 

m 

Z 2N > Tr(RR T p™) = E E H, R aA (3-5) 

ab £=i 

where R is the rectangular matrix (R\, R2, ■ ■ ■ , Rm)- The matrix product RR T is no more an 
identity matrix, but is a kind of projection operator that satisfies 

RR T =(RR T y, (3.6) 

where m = r&nk(RR T ) coincides with Tr(^RR T ^j. Let us consider a g^-dimensional matrix I that 
satisfies I = I 2 and rank(J) = m. One finds Tv(^RR T p 2N ^j > Tr(^Ip 2N ^j , where the r.h.s. coincides 
with the l.h.s. only when I = RR T . If we regard Eq.(3.5) as a variational relation for the partition 
function, the projection operator RR T in eq.(3.5) gives maximum variational partition function (or 
minimum free energy) under the constraint rank (RR T ) = m. 

We can generalize the variational relation eq.(3.5) to finite size systems with arbitrary shape. For 
example, let us consider an odd-size (L = 2N + 1) square cluster, which appears in eq.(2.13) and 
eq.(2.15). We choose the example because it is a typical system whose density matrix is asymmetric. 
The density matrix is defined as 

Pah" = E W cdef (P c N C N P d N C N P e N C N P?C N ) ab , (3.7) 

cdef 

where arrangements of the spin indices are shown in Fig. 6. We have defined p 2N+1 so that the 
matrix dimension coincides with that of p 2N in eq.(3.1). Since p 2N+l is asymmetric, the formulation 
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Fig. 6. Density matrix p 2N+1 for N = 3. Unlike p 2N , the p 2N+1 is not symmetric. 

in eq.(3.2) and eq.(3.3) should be modified; we have to consider the left eigenvalue problem 

Of P ^ = \Oj (3.8) 

independently from the right one 

P 2N+1 Q i =Q i \, (3.9) 

where Aj is the eigenvalue in the decreasing order Ai > A2 > • • • > 0, and Oi and Qj are the left 
and the right eigenvector of p 2N+l , respectively. We express the vector element of Oi and Q i as 
O ai and Q ai , respectively. This time the square matrices Q = (Qi , Q21 ■ ■ ■) an d O = (0±, O2, ■ ■ •) 
are not orthogonal by themselves, while they still satisfy the dual orthogonal relation 



1 Q).. = (O i ,Q j ) = 6 ij . (3.10) 

In other word, T is the inverse of Q. The projection operator for p 2N+1 , which corresponds to 
RR T in eq.(3.5), is then given by 

m 

(QO T ) ah = Y.QaP hc (3.H) 

where Q = (Qi,Q2, ■ ■ ■ , Q m ) and O = (Oi, O2, • • • , O m ) are rectangular matrices. The dual or- 
thogonal relation eq.(3.10) ensures that the matrix O t Q satisfies (QO T ) = QO T and Tr(QO T ) = 
m. The inequality in partition function (eq.(3.5)) is modified to 

Z 2N+1 > Tr (Qd T p 2N +^ (3.12) 
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for the odd-size system. 

We have obtained two projection operators, RR T in eq.(3.5) and QO T in eq.(3.11), where they 
have the same matrix dimension (= m). Since RR T gives the maximum variational partition 
function when it is applied to p 2N , we obtain the relation 

Tt(rR t P 2N ) > Ti(Qd T p 2N ) , (3.13) 

where the r.h.s. approaches to the l.h.s. with increasing m. Conversely, we have 

Tr(QO T p 2N+1 ) >Tt(rR t p 2N+1 ). (3.14) 

Let us compare the even-size system in Fig. 5 with the odd-size one in Fig.6. Both of them have 
the same cut (or gap) of length N = 3, but their system sizes are different. The inequalities 
eq.(3.13) and eq.(3.14) show that each system has its own optimal projection operator, that is not 
the optimal one of other systems. The ratio 

Tr(M T p 2N+1 )/Tr(Qd T p 2N+1 ) (3.15) 

depends on how the density matrix catches the boundary effect. If the system is off-critical (or 
massive), the boundary effect is expected to disappear in the thermodynamic limit N — ► oo, and 
therefore the ratio eq.(3.15) approaches to unity with increasing N. Even for a system at criticality, 
it is numerically observed that the ratio approaches to unity with increasing N. 

§4. Renormalization group algorithm 

The projection operator RR T in eq.(3.5) restricts the degree of freedom of p 2N down to m. 
Therefore the operation of the matrix R on p 2N can be regarded as a RG transformation. In this 
sense the m-dimensional diagonal matrix 

p 2N = R T p 2N k = dfcg^u* . . . ^J. (4.1) 

is the renormalized density matrix, that satisfies the variational relation 

Z 2N > Tr ~2N = ^(My") • (4-2) 

Since p 2N is written as a product of four CTMs, the RG transformation in eq.(4.2) can be naturally 
extended to CTM: 

C N = R T C N R = diag(wi,W2, . . . ,u m ) . (4.3) 

There are several ways to define the RG transformation for HRTM. A way is to use the variational 
relation for z 2N+1 in eq.(3.12). Since p 2N+1 contains as shown in eq.(3.7), the rectangular 
matrix O in eq.(3.12) transforms as T P^O. However, we don't follow this way, because 
the renormalized HRTM thus defined is not consistent with the renormalized CTM in eq.(4.2). In 
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order to avoid the problem, we define the RG transformation for HRTM using another variational 
relation 

Z 2N+1 > T ^ M T p 2N+l^ (44) 

which is derived from eq.(3.12) and eq.(3.14). This time R touches , which exists inside p 2N+1 . 
(See Eq.(3.7)) ^From eq.(4.4), we can define the renormalized HRTM as 

P? = R T P?R, (4.5) 

where the (£77) element of P^ is P^ ar) - 



Fig. 7. Extension of the renormalized CTM. Note that C N is a diagonal matrix. 



Now we explain the detail of our RG algorithm, that extends the system size recursively using 
a mapping from C N and P^ to C N+1 and P^ +1 , respectively. Remember that C N is defined 
recursively in eq.(2.8). It is possible to generalize the recursion relation to the renormalized CTM 
in eq.(4.3) 

^(S)(/3,b) = Yl W efba P ae5 P j3f5 C 55 > ( 4 - 6 ) 
efS 

where C^^p ^ is an extended CTM of linear size N + 1, and the pairs of indices (a, a) and (/3, b) 
represent the row and the column indices of the gm-dimensional matrix C N+1 . (Fig. 7) We have 
used the fact that C N is diagonal: C^p = & a f3U a . The matrix C N+1 is 'partially diagonalized' in 
the sense that it contains m-state block-spin indices a and (3. We then perform the renormalization 
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group transformation on C N+1 and decrease its matrix dimension (= qm) down to m. For this 
purpose we create the density matrix for the extended system 

, 1 



P 



2(N+1) = ( C N+1 ) (4?) 



and solve the eigenvalue problem 



?( N +VR i = R i \ i (4.8) 

in order to create a new RG transformation matrix R = (R\, R2, • • • , R m ), where the eigenvalues 
Aj are in the decreasing order. As we have done in eq.(4.3), we perform RG transformation for 
C N+1 as 

C N+1 = R T C N+1 R, (4.9) 

to obtain the new renormalized CTM C N+1 , which is an m-dimensional diagonal matrix. We also 
perform the RG transformation for HRTM at the same time. First we increase the length of HRTM 
using the relation 

d=l 

which is a generalization of eq.(2.8), and then renormalize P b N+1 by applying R to P b N+1 : 

P? +1 = R T P* +1 R. (4.11) 

In this way, we have obtained C N+1 and P N+l from C N and P N through eqs.(4.6)-(4.11). This 
is a cycle in CTMRG. Since the matrix dimension of C N in eq.(4.6) is equal to (or less than) 
qm, we can repeat the RG process up to arbitrary N. Compare to Baxter's methodic SB, the 
process of the renormalization group transformation (eq.(4.9)) is the same, but the way to extend 
CTM (eq.(4.6)) is different; the chief difference is that the renormalized HRTM appears explicitly 
in CTMRG, while it is absent in Baxter's method. 

The maximum matrix element in C and P N grows exponentially with respect to N, since free 
energy is extensive. Therefore we should take an appropriate normalization for both C and P N 
during the iteration. The simplest way is just to divide these matrices by their largest element, 
and store the normalization factor. 

Every after the iteration, the variational (or the approximate) partition function is obtained from 
the relations Z 2N = Tr^C*^ , or 

kheb 

Thermodynamic functions can be calculated from the numerical derivative of the variational free 
energy F = —ksT 'In Z 2N with respect to the temperature. Spin polarization and correlation 
functions can be expressed as a product of renormalized CTM and HRTM. For example, spin 



E W kheb Tr(p k N C N P^C N P e N C N P b N C N ) . (4.12) 
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polarization at the center of the odd-size system is calculated as 

iNfjN pNriN pN nN pNrjN 

: 'A' 



(<^is) = 7 (4.13) 



which is the renormalized expression for eq.(2.14). As shown in eq.(4.12) and eq.(4.13), the renor- 
malized system of size L = 2N + 1 has original g-state spin variables at the center. Therefore local 
quantities at the center are most easily obtained. 

The algorithm of CTMRG is closely related to the infinite system algorithm in dmrgS Okunishi 
have investigated their common thermodynamic limit© The major difference between DMRG and 
CTMRG is that the DMRG create density matrix using an eigenvector of the row-to-row transfer 
matrix, while in CTMRG the density matrix is expressed as a product of four CTMs. In other 
word, DMRG treats an infinitely long d = 2 system of width L, while CTMRG treats a square 
shaped finite size system of size L. 

We have treated the symmetric g 4 -vertex model as a reference system. It is straightforward to 
apply the CTMRG to other systems, such as the J1-J2 Ising model,0il'iil' the IRF model, and 
triangular systems. Speaking more generally, CTMRG is applicable to periodic d = 2 classical 
lattice systems that have short-range interactions and discrete spin (or site) variables. Since the 
CTM of these models are not always symmetric, care must be taken for the diagonalization of the 
density matrix; we don't explain the detail here because it is rather lengthy.!!^ 



Fig. 8. A cut of length 3 on the d — 2 lattice. The cut touches the shaded 3 by 3 square block. 
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In closing this section, we comment on the relation between the conventional block RG& and 
CTMRG. Let us consider an infinitely large square lattice that has a cut of length 3. (Fig. 9) As 
we have defined the density matrix for systems in Figs. 5-6, we can define g 3 -dimensional density 
matrix that corresponds to the cut (or gap) in Fig. 9. The density matrix naturally leads a RG 
transformation R as we have mentioned in eq.(4.2), where R transforms a set of 3 spins into an 
m-state block spin variable. By applying the block spin transformation for every side of the shaded 
square, we can map the 3 by 3 cluster into a m 4 -vertex. The CTMRG can be interpreted as a kind 
of block RG, where the block size is equal to a quadrant of the whole system. 

§5. Numerical Results 

The decorated Potts model, which we have treated as a g 4 -vertex model, has one-to-one corre- 
spondence to the Potts model on the simple square lattice; see eqs.(2.1)-(2.3). Hereafter we refer 
to the latter simply as 'Potts model'. Using the correspondence, we perform numerical calculations 
for (? 4 -vertex model and observe the scaling behavior of the q = 2, 3 Potts model at criticality 

K c = In (y/q + 1). (5.1) 

We examine the calculated data using finite size scaling and compare the evaluated critical 

exponents r/ and v with exact ones. 
The order parameter of the g-state Potts model is defined as0 

M ^ q{5l q s) _~ l . (5.2) 

At the center of square cluster, it is expected that the local order parameter M obeys the scaling 
formula 

M ~ L-^+W 2 (5.3) 

at criticality, where d is the spatial dimension (= 2), and L is the linear dimension of the system. 
We impose fixed boundary conditions in eq.(2.11), because otherwise (Si s ) is always zero. We 
evaluate (Si s ) using eq.(4.13). 

Figure 9 shows the system size dependence of the order parameter M in eq.(5.2) of the q = 2 
Potts model, i.e., the Ising model. We plot representative data when m = 4 and m = 200. 
The parameter M is almost proportional to L -1 / 8 , which is consistent with the finite size scaling 
behavior in eq.(5.2) with rj = 1/4. In fact, the least-square fitting to the calculated data in the 
range 10 < L < 1000 gives r] = 0.2504. We also perform the same scaling analysis for the q = 3 
Potts model. The obtained r\ when m = 200 is summarized in Table I. The calculated exponents 
are in accordance with the exact results. 
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Fig. 9. Log-scale plot of the order parameter M in eq.(5.1) of the q = 2 Potts model. 



Table I. Critical exponents rj of the q — 2,3 Potts Models. Least-square fittings are done in two regions, 10 < L < 
1000 and 100 < L < 1000. 



Exact 



L = 10 ~ 1000 



100 ~ 1000 



1/4 
4/15 



0.25 
= 0.26 



0.2504 
0.2652 



0.2501 
0.2654 



We calculate another critical exponent v using the finite size scaling behavior of the local energy 

E-E c ~ L l > u - d (5.4) 

at criticality, where E is the nearest neighbor spin correlation function 

E = (5(si,s i+1 )), (5.5) 

and E c is the local energy per site in the thermodynamic limit L — > oo. The value of E c is equal 

E c = \ + -L (5.6) 
2 2^q 
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where the value comes from the duality relation .0 Table II shows the calculated exponents v for 
q = 2,3 Potts models when m = 200. The calculated exponents agree with the exact ones.S'H' 

One of the recent Monte Carlo (MC) simulation up to N = 512 by Swendsen-Wang algorithm!!) 
gives the exponent v = 0.835(2) for the case q = 3© Therefore, the numerical precision of 
the CTMRG method is comparable to that of the recent MC simulation. The superiority of the 
CTMRG is that the computation time is proportional to N, while it is proportional to TV 2 in MC 
simulation. In addition, it should be noted that MC simulation requires several independent runs 
in order to collect scaling data; in CTMRG one can obtain the scaling data at once by single run. 



Table II. Critical exponents v of the q — 2,3 Potts Models. 



q Exact L = 10 ~ 1000 L = 100 ~ 1000 



2 1 1.0017 1.0006 

3 5/6 = 0.83 0.8323 0.8321 



We finally compare the numerical convergence in CTMRG to that in Baxter's variational 
methodJiBil'il) Figure 10 shows the number of iterations that are necessary for obtaining site 
energy of q = 2 potts model; we stop the iteration when the second largest eigenvalue of the 
renormalized CTM converges down to the precision 10 -8 under the condition m = 20. The thick 
line shows the iteration number by CTMRG, and the thin line shows that by Baxter's variational 
method. Throughout the whole temperature region, CTMRG exhibits better numerical conver- 
gence than Baxter's method. In particular, at the temperature shown by the cross marks, the 
calculation by Baxter's method is 'trapped' at a quasi stable point, and therefore the obtained 
results are not correct; such an instability does not occur in CTMRG. 

§6. Conclusion 

We have explained the minimum free energy principle, which is the common back ground in both 
DMRG and Baxter's CTM method. ^From the variational view point, we introduce the concept 
of DMRG into Baxter's method, and obtain a new RG method for d = 2 classical lattice models. 
Apart from the original DMRG, the numerical algorithm of our method (CTMRG) is stable even 
at the critical point. 

In order to check the efficiency of CTMRG, trial calculations are performed for q = 2, 3 Potts 
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Fig. 10. Number of numerical iterations that are necessary for obtaining site energy. The thick line shows the 
number by CTMRG, and the thin line shows that by Baxter's variational method. At the point shown by the 
marks, the free energy did not obtained correctly by Baxter's method. 



models at criticality. Calculated data are analyzed using the finite size scaling method, and it is 
confirmed that critical exponents of these models are correctly evaluated by CTMRG. 

The authors would like to express their sincere thanks to Y. Akutsu and M. Kikuchi for valuable 
discussions. T. N. thank to G. Sierra, M. A. Martin-Delgado, and S. R. White for helpful discussions 
about DMRG. K. O. is supported by JSPS Research Fellowships for Young Scientists. The present 
work is partially supported by a Grant-in-Aid from Ministry of Education, Science and Culture 
of Japan. Most of the numerical calculations were done by NEC SX-3/14R in computer center of 
Osaka University. 
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